#libraries and modules
import math as m
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats.stats as st

sin=m.sin
cos=m.cos
rad=m.radians
tan=m.tan



def infiniteslope_FS_seep(c1,gamma,z,B,ru,phi1):
    return (c1/(gamma*z)  +  (cos(rad(B))**2 -ru)*np.tan(np.radians(phi1))) /  (cos(rad(B))*sin(rad(B)))

def infiniteslope_FS_EQ(c1,gamma,z,B,phi1,kh):
    kv=0.5*kh
    c1,phi=strength_red_EQ(kh,c1,phi1)    
    return (c1/(gamma*z)  +  ((1+kv)*cos(rad(B))**2 -kh*cos(rad(B))*sin(rad(B)))*np.tan(np.radians(phi1))) /  ((1+kv)*cos(rad(B))*sin(rad(B))+kh*cos(rad(B))**2)


def strength_red_EQ(kh,c,phi):
    PGA=2*kh*9.81
    red=np.where(PGA<0.01,
                 1*PGA,
                 np.where(PGA<0.05,
                          ((0.10-0.01)/(.05-.01)) *  PGA + (0.1-0.05*((0.1-0.01)/(.05-.01))),
                          ((0.50-0.10)/(.40-.05)) *  PGA + (0.5-0.40*((0.5-0.10)/(.40-.05)))))
    
    return c*red, phi*red
        
        


def unifdist(mn,mx,n):
    return np.random.random(n)*(mx-mn)+mn


#MAIN
loop='y'
while loop=='y' or loop=='Y':

    print "FACTOR OF SAFETY ANALYSIS FOR INFINITE SLOPES USING EFFECTIVE STRESS\n"
    print "INPUT PARAMETERS\n"


    #slope gradient in degrees
    B=float(raw_input("Slope gradient (degrees): "))
    H=float(raw_input("Slope height (m): "))



    #thickness of soil in m
    repeatz=1
    while repeatz==1:
        max_z=float(raw_input("\nMaximum soil thickness (m); must be less than "+str(round(H/3.,2))+": "))
        min_z=float(raw_input("Minimum soil thickness (m): "))
        if max_z>H/3. or min_z>H/3.:
            print "\nERROR: Soil thickness values not applicable for infinite slope geometry."
        else:
            repeatz=0
                    
    
    #unit weight in KN/m3
    min_gamma=float(raw_input("\nMinimum gamma (KN/m3): "))
    max_gamma=float(raw_input("Maximum gamma (KN/m3): "))

    #friction angle in degrees
    min_phi1=float(raw_input("\nMinimum phi (degrees): "))
    max_phi1=float(raw_input("Maximum phi (degrees): "))

    #cohesion in KPa
    min_c1=float(raw_input("\nMinimum cohesion (KPa): "))
    max_c1=float(raw_input("Maximum cohesion (KPa): "))


    #Ru
    print "\nFor seepage analysis only:"
    min_ru=float(raw_input("\nMinimum Ru (0-0.7): "))
    max_ru=float(raw_input("Maximum Ru (0-0.7): "))

    #kh
    print "\nFor earthquake analysis only:"
    min_kh=float(raw_input("Minimum kh (0-0.2): "))
    max_kh=float(raw_input("Maximum kh (0-0.2): "))

    #number of iterations
    num_iter=5000#int(raw_input("\nNumber of iterations (at least 1000): "))



    #generating uniform distributions
    gamma_list=unifdist(min_gamma,max_gamma,num_iter)
    phi1_list=unifdist(min_phi1,max_phi1,num_iter)
    c1_list=unifdist(min_c1,max_c1,num_iter)
    z_list=unifdist(min_z,max_z,num_iter)
    ru_list=unifdist(min_ru,max_ru,num_iter)
    kh_list=unifdist(min_kh,max_kh,num_iter)

    

    #solving for FS
    FSseep=infiniteslope_FS_seep(c1_list,gamma_list,z_list,B,ru_list,phi1_list)

    FSEQ=infiniteslope_FS_EQ(c1_list,gamma_list,z_list,B,phi1_list,kh_list)

    FS_list=(FSseep,FSEQ)

    a=0
    print "\nRESULTS\n"
    for FS in FS_list: 

        if a==0:
            analysis='seepage, ru'
        else:
            print "\nEarthquake analysis not yet available."
            break
            analysis='earthquake, kh'

        a=a+1 
        #plotting histogram of FS
        plt.figure()
        plt.hist(FS,bins=25,normed=True, cumulative=False, histtype='bar',)
        
        stats1_0=st.percentileofscore(FS,1.0)
        stats1_2=st.percentileofscore(FS,1.2)
        plt.title('Factor of safety analysis (with '+analysis+')\n%FS<1: '+str(round(stats1_0))+'       %FS<1.2: '+str(round(stats1_2)))
        plt.xlabel("Factor of safety")
        plt.ylabel("normed frequency")
        print analysis+" analysis"
        print '%FS<1.0: ',round(stats1_0,1)
        print '%FS<1.2: ',round(stats1_2,1),
        if round(stats1_2,1)>0.05:
            print "; Failure by infinite slope mechanism is admissible."
        else:
            print "; Failure by infinite slope mechanism is NOT admissible."
        print 'mean FS: ', round(np.mean(FS),1) 


        #parameters where FS~1:
        tol=0.01   #tolerance value for FS=1
        ind_les1=np.where((FS<=1-tol))
        ind_gre1=np.where((FS>=1+tol))
        ind_eq1=np.where((FS>1-tol)*(FS<1+tol))

        fig,axes=plt.subplots(nrows=2, ncols=2,sharex=True)
        
        plt.suptitle('Factor of safety analysis (with '+analysis+')\n%FS<1: '+str(round(stats1_0))+'       %FS<1.2: '+str(round(stats1_2)))
        

        curax=axes[0,0]
        plt.sca(curax)
        curax.plot(phi1_list[ind_gre1],c1_list[ind_gre1],'gx',label='FS>1')
        curax.plot(phi1_list[ind_les1],c1_list[ind_les1],'r+',label='FS<1')        
        curax.plot(phi1_list[ind_eq1],c1_list[ind_eq1],'yo',mew=1.5,label='FS=1')
        curax.set_ylabel("cohesion, KPa")
        plt.setp(curax.get_xticklabels(),fontsize='small')
        plt.setp(curax.get_yticklabels(),fontsize='small')
                 
                 
        curax=axes[1,0]
        plt.sca(curax)
        curax.plot(phi1_list[ind_gre1],z_list[ind_gre1],'gx',label='FS>1')
        curax.plot(phi1_list[ind_les1],z_list[ind_les1],'r+',label='FS<1')        
        curax.plot(phi1_list[ind_eq1],z_list[ind_eq1],'yo',mew=1.5,label='FS=1')
        curax.set_xlabel("phi, degrees")
        curax.set_ylabel("soil thickness, m")
        plt.setp(curax.get_xticklabels(),fontsize='small')
        plt.setp(curax.get_yticklabels(),fontsize='small')
                 
                 
        

        curax=axes[0,1]
        plt.sca(curax)
        curax.plot(phi1_list[ind_gre1],gamma_list[ind_gre1],'gx',label='FS>1')
        curax.plot(phi1_list[ind_les1],gamma_list[ind_les1],'r+',label='FS<1')        
        curax.plot(phi1_list[ind_eq1],gamma_list[ind_eq1],'yo',mew=1.5,label='FS=1')
        plt.legend()
        curax.set_ylabel("gamma, KN/m3")
        plt.setp(curax.get_xticklabels(),fontsize='small')
        plt.setp(curax.get_yticklabels(),fontsize='small')

        curax=axes[1,1]
        plt.sca(curax)
        curax.plot(phi1_list[ind_gre1],ru_list[ind_gre1],'gx',label='FS>1')
        curax.plot(phi1_list[ind_les1],ru_list[ind_les1],'r+',label='FS<1')        
        curax.plot(phi1_list[ind_eq1],ru_list[ind_eq1],'yo',mew=1.5,label='FS=1')
        plt.setp(curax.get_xticklabels(),fontsize='small')
        plt.setp(curax.get_yticklabels(),fontsize='small')
        curax.set_xlabel("phi, degrees")
        curax.set_ylabel("ru")

        fig.tight_layout()
        fig.subplots_adjust(top=0.9)
        
        
        

        t=0
        if t==1:
            fig,axes=plt.subplots(nrows=5,ncols=1,sharex=True)
            axes[0].plot(FS,z_list,'.')
            axes[0].set_ylabel("Z, m", rotation='horizontal')
            plt.setp(axes[0].get_yticklabels(), visible=False)
            
            axes[1].plot(FS,gamma_list,'.')
            axes[1].set_ylabel("gamma, KN/m3", rotation='horizontal')
            plt.setp(axes[1].get_yticklabels(), visible=False)

            axes[2].plot(FS,phi1_list,'.')
            axes[2].set_ylabel("phi, degrees", rotation='horizontal')
            plt.setp(axes[2].get_yticklabels(), visible=False)

            axes[3].plot(FS,c1_list,'.')
            axes[3].set_ylabel("coh, KPa", rotation='horizontal')
            plt.setp(axes[3].get_yticklabels(), visible=False)

            axes[4].plot(FS,ru_list,'.')
            axes[4].set_ylabel("ru", rotation='horizontal')
            plt.setp(axes[4].get_yticklabels(), visible=False)

            plt.tight_layout()

        

    plt.show()

    


    try:
        loop=raw_input("Input another set of parameters? <y/n>: ")
    except:
        print "\nInput error...closing..."
        break
    
